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In this paper, we study the Hubbard model with intersite Coulomb interaction in the ionic limit 
(i.e. no kinetic energy). It is shown that this model is isomorphic to the spin-1 Ising model in 
presence of a crystal field and an external magnetic field. We show that for such models it is 
possible to find, for any dimension, a finite complete set of eigenoperators and eigenvalues of the 
Hamiltonian. Then, the hierarchy of the equations of motion closes and analytical expressions for 
the relevant Green's functions and correlation functions can be obtained. These expressions are 
formal because these functions depend on a finite set of unknown parameters, and only a set of 
exact relations among the correlation functions can be derived. In the one-dimensional case we 
show that by means of algebraic constraints it is possible to obtain extra equations which close the 
set and allow us to obtain a complete exact solution of the model. The behavior of the relevant 
physical properties for the ID system is reported. 
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I. INTRODUCTION 



i 

In a recent papeni we have shown that there is a large class of fermionic systems for which it is possible to find a 
complete set of eigenoperators and eigenvalues of the Hamiltonian. Then, the hierarchy of the equations of motion 
closes and analytical expressions for the Green's functions (GF) can be obtained. 

In this article, we apply this formulation to the extended Hubbard model, where a nearest-neighbor Coulomb 
interaction term is added to the original Hamiltonian. This model is one of the simplest models capable to de- 
scribe charge ordering in interacting electron systems, experimentally observed in a variety of systems. We will 
study the model in the ionic limit, where the kinetic energy is neglected with respect to the local and intersite 
Coulomb interactions. Among the many analytical methods used to study the extended Hubbard model we recall: 
Hartree-Fock approximation^, perturbation theory^, dynamical mean field theory^, slave boson approach 5 ^, coherent 
potential approximation^. Numerical studies by means of Quantum Monte Carlo&, Lanczos technique 9 " and exact 
diagonalizationiS have also to be recalled. 

As it will be shown in Section 2, the extended Hubbard model in the ionic limit is isomorphic to the spin-1 Ising 
model in presence of a crystal field A and an external magnetic field h. The latter model is known as the Blume-Capel 
(BC) modelii*i2iiii4. With the addition of a biquadratic interaction K it is known as the Blume-Emery-Griffiths 
t— I ' (BEG) model 1 — and has been largely applied to the study of fluid mixtures and critical phenomena. The model 
is also related to the three-component modeli&. Some exact results for the BC and BEG models are known. In 
one dimension and zero magnetic field, the spin-1 Ising model and the BEG model have been solved exactly by 
means of the transfer matrix methodii±&, and by means of the Bethe methodic. Exact solutions have also been 
obtained for a Bethe latticed and for the two-dimensional honeycomb latticed. The most common approach to the 
BC and BEG models is based on the use of mean field approximationi 5 *iii*22i22i24i2 5 ri2£i2I. However, renormalization 
. group Studies22i22i22i2ii22i2ii2i show some qualitative differences from the mean field results. Among other techniques, 
we mention temperature expansionsiSSSiS, cluster- variation method^ and numerical simulations^iMii^. The self- 
consistent Ornstein-Zernike approximation has been used to study the phase diagram of the 3D Blume-Capel model for 
O ■ spin 1— and spin 3/2-i A generalization of the BEG model by introducing a nonsymmetric exchange interaction L was 
introduced in Ref. l23l The one-dimensional case for this model was studied in Ref. liM where exact renormalization- 
group recursion relations were derived, exhibiting tricritical and critical fixed points. Also, it should be mentioned 
that the general spin-1 model can be mapped onto the spin-1/2 Ising model, under certain constrained conditions, 
which determine the corresponding subspaces of interaction parameters ( J, L, K, h, A)2£. The Ising model with spin 
1/2, 1 and with general spin has been also studied by means of the Green's function f 0lma li S x n gJ&J2£2£i£2£2£££!ii£&. 
These studies do not lead to a complete solution, but to a series of exact relations among the spin correlation functions. 
These correlation identities have been used as basis for high temperature expansions'* 5 ^, and in combination with 
the effective- field approximation 58,59 . 

The outline of the paper is as follows. In Section 2, we introduce the model for a cZ-dimensional cubic lattice. 
In Section 3, we show that it is possible to find a closed set of composite operators, which are eigenoperators of 
the Hamiltonian and close the algebra. Then, as shown in Section 4, analytical expressions of the retarded Green's 
function (GF) and correlation function (CF) can be obtained. These expressions are only formal. As the composite 
operators do not satisfy a canonical algebra, the GF and CF depend on a set of internal parameters, not calculable 
through the dynamics, and only exact relations among the correlation functions are obtained. In the framework 
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of the Green's functions formalism, extra equations must be found by fixing the representation. According to the 
scheme of the composite operator methodffiiSiiS (COM), we fix the representation by means of the algebra (algebra 
constraints). By following this scheme, in Section 5 we are able to derive for the one-dimensional case extra equations 
which close the set of relations and allow us to obtain an exact solution of the ID extended Hubbard model in the 
ionic limit. This solution is also a solution of the ID spin-1 Ising model in presence of a crystal field and an external 
magnetic field. As already mentioned, by using the GF f orma li sm 2Li&2SL5ili5ii52i§£5iiSi5£, other authors have derived 
a set of exact equations for the correlation functions of the Ising model. All of them did not succeed to find the extra 
equations necessary to close the set. There is one exception: in Ref. l53l the set of equations for the ID spin-1/2 Ising 
model for an infinite chain is closed by using ergodicity conditions for the correlation functions. However, it should 
be remarked that ergodicity breaks down for finite systems and at the critical points. In Section 6, we present some 
results for the particle density, specific heat and compressibility, both in the case of attractive and repulsive intersite 
Coulomb interaction. Details of the calculations are given in the Appendices. 

We would like to comment that the extended Hubbard model, although in the limit of localized electrons, is of 
physical interest. The results reported in Section 6 show some relevant features: (a) the behaviors of the particle 
density and of the double occupancy show the occurrence, at T = 0, of phase transitions towards charge ordered 
states (in particular, for V > 0, a checkerboard order establishes in the region < /i < 4V); (b) the specific heat 
presents a double peak structure; (c) a crossing point in the specific heat curves can be observed (it is remarkable to 
note that this crossing is observed only in the region where the checkerboard order is present and the compressibility 
vanishes); (d) in the low T region, the thermal compressibility exhibits a double peak structure, with peaks localized 
at /i = and fj, = AV. All the above mentioned results are characteristic of the Hubbard interactions and are somehow 
independent of the mobility of the electrons. Indeed, very similar results have been obtained for the complete Hubbard 
model (i.e., with finite hopping) by making use of approximations. 

II. THE MODEL 

A simple generalization of the Hubbard model is obtained by including an intersite Coulomb interaction. The 
Hamiltonian of this model is given by 

H = X>« - ^n]c^i)c(j) + Uj2n r (i)n L (i) + \ £ V V} n{i)n{j) (1) 
with the following notation. c(i) and c'(i) are annihilation and creation operators of electrons in the spinor notation 

c «=(cl(i)) ct(i) = ( 4(i) c{(i)) (2) 
and satisfy canonical anti-commutation relations: 

{c CT (i,*),4'(j>*)} = 5^5^ (3) 

{c-CMwa.*)} = {4(i,t),4,(j,t)}=o 

i stays for the lattice vector Ri and i = (i, t). The spinor notation will be used for all fermionic operators, [i is the 
chemical potential, iy denotes the transfer integral and describes hopping between different sites. n a (i) = cJ.(i)c CT (i) 
is the charge density of electrons at the site i with spin a. The strength of the local Coulomb interaction is described 
by the parameter U . n(i) is the total charge density operator 

n(i) = X)4(i)ca(i) = ct(i)c(i) (4) 

(7 

and Vlj describes the intersite Coulomb interaction. In this work, we restrict the analysis to the ionic limit (i.e., 
ty = 0). By considering only first-nearest neighboring sites, Vy = — 2dVay where d is the dimensionality of the 
system and ay is the projection operator. For a cubic lattice of lattice constant a we have 



a(k) = - V* cos(fc„a) (5) 
d ^— ' 

n— 1 



k 
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Then, the Hamiltonian (JTJ takes the form: 

H = ^[-/zn(i) + UD{i) + dVn{i)n a {i)] (6) 

i 

where we introduced the double occupancy operator 

D{i) = n,(i) ni {i) = l -n{i)[n{i)-l] (7) 
Hereafter, for a generic operator we use the following notation 

$ a (M)=5> y <J>Ci,f) (8) 



Let us consider the transformation 



It is clear that 



i(i) = l + S(i) (9) 



n(i) = ^ S(i) = -1 

n(i) = 1 S(i) = (10) 

n(i) = 2 5(i) = 1 

Under the transformation © the Hamiltonian © can be cast in the form 

H = -djY,S(i)S a (i) + S 2 (i) -h^Sii) + E (11) 



where we defined 



£ = (-ju + rf^)^ J = -dV 
h = n- 2dV -iU A = \U 



(12) 



Hamiltonian (|ll|l is just the spin-1 Ising model with first-nearest neighbor interactions in presence of a crystal field 
A and an external magnetic field h. We have the equivalence 

Hlsing = H — E (13) 

The relation between the partition functions is 

Z„=e-P E °Z Ising (14) 
then, the thermal average of any operator A assumes the same value in both models 

Wh ^ ( A )lsin 9 (15) 

According to this we can choose to study either one or the other model and obtain both solutions at once. We decided 
to put attention to the model Hamiltonian (JBJ. However, all the results can be easily extended to the Ising model 
by means of the property l|15(l and of the transformation rules © and (|12J) . In closing this Section, we note that 
the particle-hole symmetry enjoyed by the Hubbard model corresponds to the symmetry of the Ising model under 
simultaneous inversion of spin and magnetic field. 

III. COMPOSITE FIELDS AND EQUATIONS OF MOTION 

It is immediate to see that the density operator n a (i) does not depend on time 

i w n a (i) = [n a (i),H] = (16) 
at 
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and the standard methods based on the equations of motion are not applicable in terms of this operator. In order to 
use the Green's function formalism, let us introduce the composite field operators 



^ ) (0 = »?(0[n a (i)] p - 1 



(17) 



where = [1 — n(i)]c(i) and r)(i) — n{i)c(i) are Hubbard operators in the spinor notation [see (J2J]. The field 
operators tpp (i) and ^^(i) satisfy the equations of motion 



(18) 



Apparently, the equations of motion do not constitute a closed set. By taking higher-order time derivatives we 
generate a hierarchy of composite operators. However, on the basis of the anticommutation relations the following 
fundamental properties of the field [n°(«)] p can be established 



-id 



(19) 



m— 1 

where the coefficients are rational numbers which satisfy the relation 

id 



i 



(20) 



The recurrence relation 1191) is proved in Appendix A. We now define the composite operators tp^(i) and ip^(i), 
multiplet operators of rank 4d + 1 



( m 



\ 



(21) 



( ^(0 \ 

^(0 



ri(i)[n a (i)) 



(22) 



\v(i)[n a (i)] 4d J 

By means of ()18|l and of the recurrence formula 119|) , these fields are eigenoperators of the Hamiltonian © 



t 



(23) 



where £® and are the energy matrices, of rank (4<i + 1) x (4d + 1), which can be calculated by means of the 
equations of motion i|18|) and the recurrence rule l|19f) . The explicit expressions of the energy matrices are given in 



Appendix B. The eigenvalues and E^ of the energy matrices are given by 



(€) 

n 

(v) 



-H + (n-l)V 

-(j, + U + (n - 1)V 



{n = l,2- 



•(4d+l)} 



(24) 



The Hamiltonian © has been solved since we have obtained a closed set of eigenoperators and eigenvalues. Then, 
we can proceed to the calculation of observable quantities. This will be done in the next Sections by using the 
formalism of Green's functions (GF). It is worth noticing that although at the level of equations of motion the two 
fields ip(€'(i) and %p™ii) are decoupled, they are indeed coupled by means of the self-consistent equations necessary 
to determine the correlators appearing in the normalization matrix [see Section 4]. 
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IV. RETARDED AND CORRELATION FUNCTIONS 



Let us consider the retarded Green's function 



(25) 

where (■ • •) denotes the quantum-statistical average over the grand canonical ensemble and a,b = It can be 

shown that G^ ab \i,j) = SabSijG^ (U — tj). By introducing the Fourier transform 



(27T) 



(26) 



and by means of the field equations (|23[) . the retarded GF satisfies the equation 

[w-eW]GW(w)=/W 
where /( a ) is the normalization matrix, defined as 

/(°) = /{^°>(^ (a)t (*)}! 



(27) 



(28) 



Calculations of the anticommutator for a paramagnetic phase and use of the recursion rule (|19|l show that the 
normalization matrix has the following expression 



1,1 



r(«0 - 



rO) 

1,2 

r (a) r (a) 
J l,2 z l,3 



r(«) 
^1,4^ 
r(a) 

J l,4d+1 



j(a) \ 
r(«) 



r(«) A a ) 
1 l,4d 1 l,4d+l 



r(«) 



(«) 



X 2,4<M-1 

r(a) r(a) 
'4d-l,4d+l J 4ci,4ci+1 



\ -n,4d+l ''2,4d+ 



r(») 



r(«) 



(29) 



1 4d,4d+l J 4rf+l,4rf+l 



where 



4d 



Ip,ld+1 ~ ^m d+1 ^p-l,m+l (P — 2, 3, • • • , Ad+ 1) 



(30) 



m— 1 



We see that we need to calculate only the 4d+ 1 elements l[ a p (p = 1, 2, • • • , 4d+X). These elements have the following 
expressions 



with the definitions 



jg = k^p-V - At*" 1 ) 
/^(k)-A(p- 1 ) 



k (p) = ([?i Q (i)]P) 

AW = \ (n(i)[n a (i)]P) 



The solution of Eq. J27J is 



4ci+l 



(a,n) 



where the spectral functions a^°'™' > are calculated by means of the formula^ 



(a,n) _ n \a) ^ roWl-l T\ a ) 
8 



»(«)• 



•>( a )l-l 7-( a ) 



(31) 



(32) 



(33) 



(34) 
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where is the {Ad + 1) x {Ad + 1) matrix whose columns are the eigenvectors of the matrix e^ a \ Calculations of 
the matrices QS^ and Qy 1 ' are reported in Appendix B. It is worth noting that we have fiK> = SlW. 

The spectral density matrices cr^"'™^ are calculated in Appendix C. They satisfy the sum rule 



4d+l 



This is a particular case of the general sum rule 



4d+l 

J2 [E^fa^ = (36) 

n=l 

where M^ a ' p ^ are the spectral moments defined as 

M {a ^ = ^{(id/di)V a) W,^ (a)t W}) (37) 

The fact that the sum rule 1)36(1 is satisfied at all orders in p, is a consequence of the theorem proved in Ref. [see 
also pag. 572 in Ref. Il^- The correlation functions can be immediately calculate from (|33|l by means of the spectral 
theorem and are given by 

C (o) (*< -tj) = (^(MiW^Mi)) = 7^) / due-^-^C^iu) (38) 



with 



4d+l , , , . 

C<°> (w) = 2^ £ if^U" ~ ^ (39) 



Equations (|33fl and (|39|) are an exact solution of the model Hamiltonian © . One is able to obtain an exact solution 
as the composite operators ipp {i) and ipp V \i) constitute a closed set of eigenoperators of the Hamiltonian. However, 



as stressed in Ref. I6CI the knowledge of the GF is not fully achieved yet. The algebra of the fields ip^{i) and ^^{i) 
is not canonical: as a consequence, the normalization matrix i n the equation l|27|l contains some unknown static 
correlation functions, correlators, [see Eqs. I|31l) - (|32l) ]. that have to be self-consistently calculated. According to 
the scheme of calculations proposed by the composite operator method£2i£ii§£, one way of calculating the unknown 
correlators is by specifying the representation where the GF are realized. The knowledge of the Hamiltonian and of 
the operatorial algebra is not sufficient to completely determine the GF. The GF refer to a specific representation 
(i.e., to a specific choice of the Hilbert space) and this information must be supplied to the equations of motion that 
alone are not sufficient to completely determine the GF. The procedure is the following. We set up some requirements 
on the representation and determine the correlators in order that these conditions be satisfied. From the algebra it is 
possible to derive several relations among the operators. We will call algebra constraints (AC) all possible relations 
among the operators dictated by the algebra. This set of relations valid at microscopic level must be satisfied also at 
macroscopic level, when expectations values are considered. Use of these considerations leads to some self-consistent 
equations which will be used to fix the unknown correlator appearing in the normalization matrix. An immediate set 
of rules is given by the equation 

(^ (a) (# (a)t «) = (2^) / duC^{u) (40) 

— oo 

where the l.h.s. is fixed by the AC and the boundary conditions compatible with the phase under investigation, while 
in the r.h.s. the correlation function C^ a \uj) is computed by means of the equations of motion [cfr. Eq. (|39j . 

Another important set of AC can be derivedi*S4 by observing that there exist some operators, O, which project out 
of the Hamiltonian a reduced part 

OH = OH a (41) 

When Hq and Hi = H—Hq commute, the quantum statistical average of the operator O over the complete Hamiltonian 
H must coincide with the average over the reduced Hamiltonian Hq 

Tr{Oe-P H } = Tr{Oe- pHa } (42) 
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Another important relation is the requirement of time translational invariance which leads to the condition that 
the spectral moments, defined by Eq. I|37|) . must satisfy the following relation 

Mf»M(k) = [M(,t' p) (k)]* (43) 

It can be shown that if (|43|l is violated, then states with a negative norm appear in the Hilbert space. Of course the 
above rules are not exhaustive and more conditions might be needed. 

According to the calculations given in appendix C, the GF and the correlation functions depend on the following 
parameters: external parameters (fi,T,V,U), internal parameters kW , , • • • , and AW , A< 2 > , • • • , A( 4d > . By 
means of the algebraic relations 

^H = !" nf (44) 

and by making use of the AC (|40|) , we obtain the following Ad + 1 self-consistent equations 

Cg + = k^-V - A (fe - X) (fc = l,2,---4d+l) (45) 

where, recalling and (p^ 

n=l (46) 

T^ = l+tanh(£i) 

To determine the 8<i parameters we need other 4d — 1 equations. In order to obtain a complete solution of the model, 
we must calculate these parameters. This will be done in the next Section for the one-dimensional case. 

It is worth mentioning that the formulation given in this Section can be easily extended to multipoint correlation 
functions, as (n(i)[n a (i)] f 'n(/i)n(Z 2 ) • • ■n(l s )}. Let us define the retarded Green's function 

G^Ht-t') = ^[V> ( °)(i,t)V (a)t ( i .*')*]) (47) 

where $ = ${n(j)} is any function of the n(J) with j i. All the equations derived above remain valid by means of 
the substitutions 

j(a) > j(o,*) _ ({^(«)(i),^(*)t(i)}$) 

K (p) — y k (p,*) = ([ n Q (i)]P$) (48) 
A(p) — y \(p,*) = i ( n (i)[n a {i)]P<S>) 

For each choice of the function &{n(j)}, it is necessary to determine the new set of parameters k^ p '^ and \( p >®\ For 
example see Ref.0, where the correlation function (n(i)[n a (i)] p n(j)) has been calculated for the ID spin-1/2 Ising 
model. 



V. SELF-CONSISTENT EQUATIONS FOR ONE-DIMENSIONAL SYSTEMS 

Until now the analysis has been carried on in complete generality for a cubic lattice of d dimensions. We now 
consider one-dimensional systems, and in particular we will study an infinite chain in the homogeneous phase. As 
shown in previous Section, the set of self-consistent equations (|45(l are not sufficient to determine all the 8 internal 
parameters. The remaining three equations can be derived by algebraic considerations on the basis of the requirement 
(|42|1 . We start from the algebraic relations 



e t «n(i) = , 



which imply that 



(50) 
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where 

H = H -2Vn(i)n a (i) (51) 
By means of the fact that H commutes with Hj = H — H , the relation (|50|) leads to 

£}{i)e- 0H = ^(i)e-? H ° (52) 
Then, by means of the requirement l|42l) . the correlation function = (C(*)^(*)[ nQ (*)] fc_1 ) can be expressed as 

r (A) rt«,o) 

r tt) - r (e,o) ^ 

U U °11 

where 

Cf 0) = <f(i)f t (i)[n a (i)]*- 1 > (54) 
and (•• -) Q denotes the thermal average with respect to iJo- In order to calculate c[^°\ let us define the retarded GF 



G§°\t-t') = (R[ttiMht')]K(i)] k - 1 ) 
G ( ?f{t-t>) = (R[r,(i,t)ri(i,t')][n^{)] k -^ 



By means of the equations of motion 



[r}(i),H } = -([M-U)r}(i) 



we have for an homogeneous phase 



2{uj + n - U+iS) 



Recalling the relation between retarded and correlation functions we have 



yiL 0) 2([»°j t -') -(n( i )K( i )] t - 1 ) 
2(1 + e^) 



2(1 + ePb-V)) 



Recalling the algebraic relations 



we obtain from Ij59(l and i jfiOjl 



where 



r) a r)t =n a - n^n l 



(1 + 2e^ + e 2 ^e~^ u ) 



(55) 



(56) 



{t 0) 2 ( n°(i) *- 1 ) n - (n (t K (» ) n 

G 1' (w) = — — — ^ — V V V U ^ (57) 

lfe 2(w + /u+i<5) 



G (f ) ( w) = v;^l- wi. ; /» (58) 



C«,o, = M wj /o n l wj /o (59) 



(61) 



(nWK(i)] fc - 1 >o= J B 1 (K(i)] fe - 1 > 

(i}( l )K( 8 )] fc - 1 ) = 5 2 (K( ! )f- 1 ) n [bZ) 



Bi = <n(*)) = ;q ; ^ ; ^r^flTn ( 63 ) 
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By substituting (^3) into ((S3) and (fHDfl 

rM>°) - n _ r, 4- n„Wr««Mi*-i\ 

(65) 



Cjf ) = i(i3 1 -2S 2 )<K(z)] fe - 1 ) 



By substituting the first equation of (|65[) into (|53[) we obtain 

C^ = c[f([n a (i)] k - 1 ) (66) 

Now, we observe^ that H describes a system where the original lattice is divided in two disconnected sublattices 
(the chains to the left and right of the site i) . Then, in Hq representation, the correlation functions which relates sites 
belonging to different sublattices can be decoupled: 

(a(j)b(m)) = (a(j)) (b(m)) (67) 

for i and j belonging to different sublattices. By using this property and the algebraic relation (IA4I) we have 

[n a (i)] 2 ) = \X 1 +X 2 + \Xl 



[n a {i)f) o = §*i + §* 2 + §*i* 2 + \X\ (68) 
\n a {i)f) o = |Xi + |*2 + |*i*2 + |*i 2 + |*| 



where we defined 



*i = (n a (i)) 
* 2 = (D a (i)) 



(69) 



Then, we obtain the self-consistent equations 



/~r(£) r 1 V" l 3 y i 3 ~v ~v i 3 v2 



U M 1 2 v ^ 1 2^ 11 ^ 1 - 2 1 4 
_ ,nrt5; rly- i 7 y i 9 y y i 7 y 2 , 3 y2 
°15 — °11 ^8 Al + 4 A 2 + 2 AlA 2 + jAj + 2 A 2 



(?) _ /nr(£) ri y , 7 y , 9 y v , 7 y2 , 3 y2l ( 70 ) 



which relate the correlation functions c[f\ c[f to c[f , c[ 2 , C13 5 when we observe that, by means of (JHEJl, the two 
parameters *i and X 2 are expressed in terms of the correlation functions c[f , C}| , C13 as 

(0 



-o(S) 1 r (t) 1 ^(?)2 

y _ °13 _ ^ °12 _ 1 12 / 79 \ 

2 r «) 2 2 r ( « } 2 1 ' 

O n O n 



We need another equation. To this purpose, we start from the algebraic relations 

DP{i) = D(i) 



D(i)nP(i) = 2PD(i) P ~ 1 ( 73 ) 



From here, with some effort, we can derive the following relations 

4 



D(i)e-fi H = + (2/p + 9 P )[n a (me- 0Ho (74) 

P =i 



^(z)e-' 3H = ^( i ){l + ^[/ p n( l )+. 9pJ D( l )]K(z)F}e^ (75) 

P =i 
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where 



h = 2A X - fAf - \A\ + \A X A 2 
h = f A i + T A l - 2 r A - 1 A 2 

h = -f A \ 2A\ + fA x A 2 (76) 
h = \A\ + \A\ \A X A, 
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-AA X + 2A 2 + f A\ - fA 2 - \A\ - fA x A 2 + U 2 A 4 



52 = -f A\ + fAl + %AZ + fA,A 2 - f A 2 A A 
"M? - fA 2 - 2AI - ^A,A 2 + f 

'16 42 , 4 A 2 , 2 A 2 , 16 >i 4 8 



3 1 3 2 1 6 4' 3 "i"-! 3 -"-^"i ( 77 \ 
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.94 = -f i4f + |A| + \A\ + fA 1 A 2 - IA 2 A 4 
and 

A p = e~ pf3V - 1 (78) 
By taking the expectation value and by using the relations and ( I66fl . we obtain from l|74|) and (|75|l 

<2J(i)) = (1-B 2 +B) {C $ + S ^ + 3p)C^ + i} (79) 

\ > p=i 



(£ a (0> = fr j L p J (^W)o + E ( B ^ + B ^p) (D a (i)ln a (i)] p }o\ (80) 



The translational invariance requires that (D a (i)) — (D(i)). Then, from l|79(l and H8()(l we obtain the equation 

4 £,(?) 4 

B 2 + B 2 J2 (2/p + ff^TT = X2 + E ( Bl ^ + B2 -9p) < ^(OKWF >o (81) 



C 



(0 

11 p=l 



To calculate the correlation functions (D a (i)[n Q (j)] p ) , we observe the following algebraic relation which can be 
derived by means of l|A4|l and (|73|) 

D a {i)n a (i) = D a (i) + \n a {i) - \[n a {i)] 2 + ±K(i)] 3 
i^)K(z)] 2 = + ±n»(i) - ±[n««] 2 - |K(z)] 3 + iK(z)] 4 

£> a (i)[n a (i)] 3 = L» Q (z) - \n a {i) + \[n a {i)] 2 - f|K(i)] 3 + \[n a {i)Y [ ' 

D a (i)[n a (i)] A = D a (i) - %n a (i) + f K(»)] 2 - 2f [n a (i)} 3 + ff[n a (z)] 4 

By taking the expectation value of (|82[1 with respect to Ho and by using the property (|67|l . we can express the 
correlation functions (D a (i)[n a (i)] p ) as 

(£>«(iK(i)) = X 2 + iX 1 X 2 

{D a (*)[n a (*)] 3 ) o = X 2 + §M 2 + jX 2 ^ 

<^( l )K(i)] 4 ) = ^2+ §x 1 x 2 + fx| 

Then, we can we can put equation (|81J) under the form 

6 + b x X x + b 2 X 2 + b 3 X!X 2 + b 4 X 2 + b 5 X 2 = (84) 

where 

bo = B 2 

bi = n- s 1 + \(r 2 - s 2 ) + i(r 3 - s 3 ) + |(r 4 - s 4 ) 

&2 = -1 + r 2 + §r 3 + |r 4 - si - 2s 2 - |s 3 - ^-s 4 , SC;A 

^ — 3„ _i_ 9^ 1 „ 5 „ 31 „ 137 „ 
03 — 2 r 3 + 2 ^4 - 2 S 1 — 4 S 2 8" S3 16" S4 

bi = |(r 2 - s 2 ) + |(r 3 - s 3 ) + |(r 4 - s 4 ) 

8 
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FIG. 1: The particle density n is plotted as a function of: (a) the chemical potential at various temperatures for V = —1; (b) 
the temperature at various values of chemical potential for V = — 1. 



r p = (B 1 + 2B 2 )f p + 2B 2 g p 
■s p = Bif p + B 2 g p 

Recalling [see l|71fl and l|72[l ] that the parameters X\ and X 2 are expressed in terms of the correlation functions c[\ , 

Gil , C{f , Eq. gives the needed third equation 

Summarizing, we have 8 self-consistent equations (|45|l , (|70|l and l|84|) which will determine the 8 internal parameters 
k,( 3 \ k( 4 \ \( 2 \ \^ , \^ in terms of the external parameters fi, T, U and V. The set of 5 equations in 
(|45|l is a system of linear equations. This system can be analytically solved with respect to 5 parameters and we are 
left with three parameters, which are determined by the non-linear equations H7U|) and H84|) . Once these parameters 
are known, we can calculate the correlation functions and all the properties of the system. 



VI. RESULTS FOR THE ONE-DIMENSIONAL CASE 



We now present some results for the case U = 0. This situation corresponds to the ionic Hubbard model without 
local Coulomb interaction, and to the pure spin-1 Ising model, without the crystal field. In one dimension and zero 
magnetic field, thespin-1 Ising model and the BEG model have been solved exactly by means of the transfer matrix 
method (Refs. Il7ll8l45j) . We recall that the case of zero magnetic field corresponds to the case of half filling in the 
Hubbard model. The presence of magnetic field has been treated in Ref. 0, only for a ferromagnetic coupling and 
studied only in connection to the existence of critical and tricritical fixed points. The general case of U different from 
zero will be considered elsewhere. We study the behavior of the system as a function of the parameters /i and T. We 
take \V\ = 1: all energies are measured in units of \V\. 

At first, we consider the case of an attractive intersite Coulomb potential (i.e., V < 0). This situations corresponds 
to J positive (i.e., ferromagnetic coupling). In Fig.n(a)j we show the particle density n =< n(i) >— as a function 
of the chemical potential \x. In terms of the Ising model, this figure should be read as the magnetization {S(i)) versus 
the magnetic field h. By increasing /i, the particle density increases and varies between and 2. At /i = 2V we have 
ri = 1, in agreement with the particle-hole symmetry. By decreasing the temperature, at /i — 2V the system tends to 
become unstable against a charge ordered state (ferromagnetic order in the Ising model): the particle density jumps 
from to 2. This is also seen in Fig. ^ (b), where the particle density is plotted versus the temperature for various 
values of the chemical potential. For /_i < 2V we have liniT^o n — 0, while for (i > 2V we have lim^^o n = 2. At zero 
temperature there is a phase transition at /j = 2V from a state with no particle to a fully occupied state where the 
charge assumes the maximum value. 

The double occupancy D can be computed by means of the expression 

D = \-Cll (87) 

The behavior of D is shown in Figs. |21 (a) and [21 (b), where D is given as a function of the chemical potential and 
temperature, respectively. By increasing /i, the double occupancy increases and varies between and 1. For /i < 2V 
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FIG. 2: The double occupancy D is plotted as function of: (a) the chemical potential at various temperatures for V = — 1; (b) 
the temperature at various values of chemical potential for V = — 1. 




FIG. 3: The specific heat C is plotted as a function of the temperature at various values of chemical potential for V = — 1. 



we have liniT— >o n = 0, while for // > 2V we have liiriT->o D = 1. At zero temperature there is a phase transition at 
/i = 2V from a state where all the sites are empty to a state where all the sites are doubly occupied. The behavior of 
the parameters and as functions of /i is similar to that exhibited by n; for T = these parameters at fi = 2V 
jump from to their ergodic value. 
The specific heat is given by 

c§ 

where the internal energy E can be calculated by means of the expression 

E = -fin + UD + 2VX {1) (89) 

The specific heat satisfies the property C(/x) = C(2V — fi). Therefore, we can limit the analysis to the region 
2V < n < oo (or — oo < fi < 2V ). As shown in Fig. [31 the specific heat increases by increasing T up to a certain 
temperature, then decreases and goes to zero in the limit T — > oo. Near the transition point /i = 2V , the peak is 
sharper and is situated in the low-temperature region. By moving away from // = 2V, the peak becomes broader and 
moves to high temperatures. 

The thermal compressibility n T is given by 

t 1 dn . . 

K -^T, (90) 
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FIG. 4: The derivative of the particle density with respect to the chemical potential dnjd\i is plotted as a function of the 
chemical potential at various temperatures for V = — 1. 




FIG. 5: The particle density n is plotted as a function of: (a) the chemical potential at various temperatures for V = 1; (b) 
the temperature at various values of the chemical potential for V = 1. 



In Fig. E| ^ is plotted versus the chemical potential for various values of temperature. A peak is observed at the 
transition point /i = 2V. By decreasing T, the height of the peak increases and the compressibility tends to diverge 
in the limit T — ► 0. As a function of the temperature (k t ) t _ l= 2v exponentially diverges at low temperatures and 
decreases as ^ in the limit of large T. In terms of the Ising model, Fig. 21 should be read as the spin susceptibility 
versus the magnetic field. 

Next, we consider the case of a repulsive intersite Coulomb interaction (i.e., V > 0). This case corresponds to J 
negative (i.e., antiferromagnetic coupling). In Fig. [51(a) we show the particle density n as a function of the chemical 
potential fj,. By increasing fj,, the particle density increases from zero, reaches the value 1 at fj, = 2V, and tends to 
2 for larger values of the chemical potential. When the temperature decreases some instabilities of the homogeneous 
phase appear. In the limit T — > 0, two singularities manifest: one at /i = 0, where n jumps from to 1, the other 
at fi = 4V, where n jumps from 1 to 2. In the region < fi < 4V, n exhibits a plateau centered at /i = 2V. This 
behavior is also seen in Fig. (b), where the particle density n is given as a function of the temperature. At T = 
we have two phase transitions. At /i = the system passes from a state with no charge to a state where the charge is 
distributed in a checkerboard structure. At fj, = 4V there is a second phase transition where the system passes from 
the checkerboard structure to a state where the charge is uniformly distributed. The checkerboard structure is clearly 
seen from the behavior of the parameters and X^ p \ as shown in Figs. El While the parameters have the same 
behavior as n, with two singularities at [i = and at /i = 4V, the parameters exhibit only one singularity at 
/i = 4V. The reason of this difference is related to the fact that (tW are correlation functions between the site i and 
second-nearest neighboring sites, while A^ p ' mainly relate the site i to first-nearest neighboring sites. 

The double occupancy as a function of the chemical potential is shown in Fig. [7| (a); by increasing fx, D increases 



14 



20 



15 - 



10 




20 



15 



10 




10 



(b) 




A (4) 


r 

/ 




/ 




/ 


A <3) 








A <2) 





10 



A' 



FIG. 6: The internal parameters (a) k (2) , k (3) , k (4) and (b) A (1) , A (2) , A (3) , A (4) are plotted as functions of the chemical potential 
at various temperatures for V — 1. 



D 




FIG. 7: The double occupancy D is plotted as function of: (a) the chemical potential at various temperatures for V = 1; (b) 
the temperature at various values of chemical potential for V = 1. 



from zero and tends to 1 for large values of the chemical potential. At T = 0, as also seen in Fig. (b), D has a 
discontinuity at /j, = 0, where jumps from zero to 1/2, and another discontinuity at fi — AV where jumps from 1/2 
to 1. In the region < /i < AV D has the constant value of 1/2. Again, this show the transition to a checkerboard 
structure of the charge. 

To study the specific heat, let us distinguish the two regions 2V < fi < AV and fi > AV. In the first region, see Fig- El 
(a), the specific heat increases with temperature, exhibits a peak at a certain temperature T = T±, then decreases. 
When approaches the critical value \i c = AV, see Fig.|5](b), the specific heat develops a double peak structure with 
a broad peak at higher temperature than T±. The latter temperature decreases with fj, and tends to zero for fi —> [i c . 
It is characteristic of this region the fact that all the specific curves cross at the same temperature, independently on 
the value of \i. The crossing temperature is T* w 1.1. The fact that there is a crossing point in the specific heat curves 
versus T, when plotted for different values of some thermodynamics quantities, has been observed in a large variety 
of systemsSSiSSiSl. In the second region, see Fig. [HI (c) , C exhibits a low-temperature peak at a certain temperature 
T = T2. T2 tends to zero for [i — > /i c , increases by increasing fi. Again, close to [i c there is a double-peak structure, 
which disappears when /i moves away from /i c . It is characteristic of this region that no crossing point is observed. For 
large value of T, the two specific heats (i.e. attractive and repulsive V) tend to coincide. This is because the system 
is in a homogeneous phase, where the thermal energy predominates over the Coulomb interaction. It is interesting 
to observe that the crossing point is observed only for V > and in the region < fi < AV, where the checkerboard 
order is observed. 
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FIG. 8: The specific heat C is plotted as a function of the temperature at various values of chemical potential for V = 1: (a) 
in the region 2V < n < AV; (b) near the transition point n = AV; (c) in the region /i > AV . 




FIG. 9: The derivative of the particle density with respect to the chemical potential dn/d/j, is plotted as a function of the 
chemical potential at various temperatures for V = 1. 



The thermal compressibility k t is studied in Fig. [5J where is given as a function of the chemical potential for 
various values of temperature. When T is lowered a double peak structure develops, with peaks localized at /.t = 
and /i = AV. The heights of the peaks increases by decreasing T and tends to diverge in the limit T — > 0. It is 
worth noticing that for low temperature the compressibility is very small (zero in the limit T — > 0) in the wide region 
< [i < AV, where the phase with checkerboard order of the charge is observed. Similar results have been obtained 
in Ref. l68i where the t — V model has been studied within a cluster approximation. 



VII. CONCLUSIONS 



The Hubbard model with intersite Coulomb interaction has been studied in the ionic limit (i.e., no kinetic energy). 
This model is isomorphic to the spin-1 Ising model in presence of a crystal field and an external magnetic field. A 
finite complete set of eigenoperators and eigenvalues of the Hamiltonian has been found for arbitrary dimensions. This 
knowledge allows us to determine analytical expressions of the local Green's functions and the correlation functions. 
As the eigenoperators do not satisfy a canonical algebra, the GF and the CF depend on a set of unknown parameters, 
not calculable by means of the dynamics. By using appropriate boundary conditions and algebraic relations we 
have determined these parameters for the case of an infinite homogeneous chain. Some results for the case U = 
(i.e., no local Coulomb interaction / no crystal field) have been given. The system exhibits a different behavior 
according to the sign of V, the intersite Coulomb interaction, or the sign of J, the exchange interaction. For V < 
(J > 0, ferromagnetic coupling) at // = TV (h = 0) the system exhibits a phase transition to a charge ordered state 
(ferromagnetic phase for the Ising model) at zero temperature. For positive V (J < 0, antiferromagnetic coupling), 
the system exhibits instabilities at \i = (h = —2 \ J\) and at \i = iV (h = 2 \ J\). In the limit of zero temperature a 
phase transition to a state with a checkerboard order of the charge (antiferromagnetic phase for the Ising model) is 
observed at /i = 0. This order persists up to /i = AV, where a second transition to an homogeneous charge order is 
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observed. A crossing point in the specific heat curves is observed only for V > and in the region < \i < AV, where 
the checkerboard order is observed. In the entire region < fi < 4Y, the compressibility vanishes at low temperatures. 
Further study, in particular for the case of finite U, will be presented elsewhere. 



APPENDIX A: ALGEBRAIC RELATIONS 



Let us start by observing that because of the basic anticommutating rules J2J the number n(i) = c'(i)c(i) and the 
double occupancy D(i) = n(i)[n(i) — l]/2 operators satisfy for p > 1 the following algebra 

n P(i) = n (i) + a p D(i) 

DP(i) = D(i) a p = 2 p -2 (Al) 

n p (i)D(i) = 2PD(i) 

From this algebra several and important relations can be derived. Firstly, let us consider the operator 

2d 



"(*') = i£>(im) (A2) 



n 

2d 

m— 1 



where i m are the first neighbors of the site i. A basic relation can be derived for the operator [n a (i)] p , with p > 1. 
In this and in the following Appendices we shall present results for the one-dimensional system£&. We start from the 
equation 



m=0 



p 

Because of the algebraic relations l|Al|l we obtain 

i 

2." 



K(i)] p = 4 E f « ) ^i) P ~"M* 2 r (A3) 



K(i)] p = b$Z m (A4) 



where the operators Z m are defined as 



Zi = 2n a (i) 

Z 2 = 2D a (i) + n(h)n(i 2 ) 
Z 3 = D(h)n{i 2 )+D{i 2 )n{i x ) 
Zi = D(h)D(i 2 ) 



(%>) 

and the coefficients bm have the expressions 

b[ p) = i 



4 P) = E ( £ )=2(-l + 2^) 



m— 1 



in 



m—l 



*4° = X? ( £ ) a p - m a m = 4(-l + 3 ■ 2 P ~ X - 3^ + A p ^) 



m— 1 

By solving the system (|A4|) with respect to the variables Z m> we can obtain from (|A4|) the recurrence rule 



where the coefficients A^m are defined as 





i 

2P- 




- ±6 (p) 




1 

21'- 






A 3 - 


1 

2 p- 






4 P) = 


1 

2p- 


t 3 4 





+ 3°3 4°4 J 



(A5) 



4 P) = I? f ^ ) Op-m = 3(1 -2 p + 3P- 1 ) (A6) 



Ktf^^K^r (A7) 



(A8) 
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We note that 



E A % ] = 1 



(A9) 



m— 1 



In Table 1 we give some values of the coefficients A 
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APPENDIX B: THE ENERGIES MATRICES 



The energy matrices and can be calculated by means of the equations of motion l|19f) and ()2()ll and the 
recurrence rule (|21(l . The explicit expressions are 



,(0 







2dV 



\ 2dVA\ 



(4d+l) 



2dVA 



-I* 

(4d+l) 
4d-l 
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2dV 

-V + 2dVA[ i d d+1) 



(Bl) 



Av) _ 








2dy 







U-fx 



2dVA 



(4d+l) 



4d-l 






2e2V 



U -fi + 2dVA { . 4d+1) 



where A 



V 2dVA[ id+1) 
^ are the coefficients appearing in (|19fl . In particular, for one dimension 
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The matrices fiw and have the expressions 



/i 
o 




Vo 



2 4 
2 3 
2 2 
2 
1 



(2/3) 4 
(2/3) 3 
(2/3) 2 
(2/3) 
1 



(1/2) 4 \ 
(1/2) 3 
(1/2) 2 
(1/2) 
1 / 



(B5) 



APPENDIX C: THE SPECTRAL DENSITY MATRICES 



By means of the formulas (|34() and recalling the expressions [see Appendix B] for the energy matrices, and the 
expressions ()31(l -()32 |l of the normalization matrix, we can easily calculate the spectral matrices. Furthermore, we 
note that the matrices OS® and flw are equal. Then, the spectral density matrices have similar form in terms of the 
matrix elements of the normalization matrices. Because of the recurrence relation, we need to calculate only the first 
row of the matrices. Calculations show that the matrices cr ( - a '™- ) have the following form 

(7 («,n) =S (a) r (n) n = 1, 2. • • • , 4rf + 1 (CI) 

where si a ' are functions of the elements (m = 1, 2. • • • , Ad + 1) and are matrices of rank (Ad + 1) x {Ad + 1). 
For the one-dimensional case 



Y,[ a > = l(Ql[ a > - 251$ + 351$ - 20l[ a > + Al[ a >) 

e£°> = + 19/$ - + 4ljj (C2) 

v( a ) _ 4/or(o) 7f(») j_ 7 rW or( a h 
2^4 — 3^1,2 _ '-'1,3 + U l,i ~ Z1 l.5> 

= |(-3jg + 11/g - 12/^ + 4lg) 



-doooo) 

r[ 2) m = (l 2- 1 2- 2 2- 3 2- 4 ) 

It! = ( 1 1 1 1 1 ) (C3) 

r^ = (l (2/3)- 1 (2/3)- 2 (2/3)- 3 (2/3)- 4 ) 

ri 5 ^ = (l (1/2)- 1 (1/2)- 2 (1/2)- 3 (1/2)- 4 ) 
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